x = {linspace(0,max(Sf),100),linspace(0,max(Nf),100)};
[SS,NN] = ndgrid(x{1},x{2}); 
y = alpha0 + alpha1./(Sbar-alpha2./(alpha3+NN)-SS);
surf(x{1},x{2},y)
